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We present and discuss three discontinuous Galerkin (dG) discretizations for the 
anisotropic heat conduction equation on non-aligned cylindrical grids. Our non- 
aligned scheme relies on a self-adjoint local dG (LDG) discretization of the elliptic 
operator. It conserves the energy exactly and converges with arbitrary order. The 
pollution by numerical perpendicular heat fluxes decreases with superconvergence 
rates. We compare this scheme with aligned schemes that are based on the flux- 
coordinate independent approach for the discretization of parallel derivatives. Here, 
the dG method provides the necessary interpolation. The hrst aligned discretization 
can be used in an explicit time-integrator. However, the scheme violates conserva¬ 
tion of energy and shows up stagnating convergence rates for very high resolutions. 
We overcome this partly by using the adjoint of the parallel derivative operator to 
construct a second self-adjoint aligned scheme. This scheme preserves energy, but 
reveals unphysical oscillations in the numerical tests, which result in a decreased or¬ 
der of convergence. Both aligned schemes exhibit low numerical heat fluxes into the 
perpendicular direction and are superior for flute-modes with hnite parallel gradients. 
We build our argumentation on various numerical experiments on all three schemes 
for a general axisymmetric magnetic held, which is closed by a comparison to the 
aligned hnite diherence (FD) schemes of References'’^. 
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I. INTRODUCTION 


Strong magnetic fields are known to separate the scales of turbulent motion in thermonu¬ 
clear fusion or astrophysical plasmas. Along and perpendicular to magnetic held lines the 
mean free paths are of the order of the connection length and of the gyroradius. This in¬ 
troduces a strong anisotropy into the magnetized plasma, resulting in viscosity and heat 
conduction coefficients, which are several magnitudes higher in the parallel direction than 
into the perpendicular direction. 

In order to formulate the problem we split the gradient into parallel and perpendicular parts 
via 


V=;hV||+Vx, (1) 

where Vy = 6 ■ V, and 6(x) = B/R is the unit vector in the magnetic held direction. With 
this dehnition we can split the dihusion operator into 

A/ = VV = V • (hh ■ V + Vx)/ =; All/ + A^f. (2) 

In the following we mean to study the exemplary model of the anisotropic heat conduction 
equation 

= x\\AiiT + x±A±T, (3) 

with constant parallel and perpendicular heat conduction coefficients yy and x±- Since in 
magnetized plasmas the parallel heat conduction is orders of magnitudes higher than the 
perpendicular heat conduction, we assume that x± is zero. The energy E = dV T is 
invariant over time if the surface integral fg^dA ■ bVyT = 0 on the boundary dQ of the 
domain 

The anisotropy imposes stringent requirements on the numerical method, which is com¬ 
monly solved by aligning the coordinates with the magnetic held^. However, aligned flux 
coordinates exhibit nonorthogonal and anisotropic meshes while not resolving singular points 
(e.g.: X- and 0-point). Especially near the last closed flux surface of fusion plasmas those 
points result into severe numerical problems, which could only partly overcome by addi¬ 
tional coordinate transformations^’^. Here, the shifted metric transformation^ remedies the 
nonorthogonality and the conformal tokamak coordinates® produce isotropic meshes for a 
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proper treatment of shaped nonsingular flux surfaces. 

The contrary approach is to misalign the numerical grid to the magnetic held. Conse¬ 
quently one has to tackle numerical issues, which are most signihcantly convergence loss ® 
and numerical perpendicular heat hux that exceeds the true perpendicular heat transport^. 
Hence, proper and more sophisticated numerical discretization techniques are unavoidable. 
In general we can distinguish between non-aligned and aligned discretizations. The former 
methods have proven successful for hnite differences^ and higher-order hnite elements®’® even 
for huge levels of anisotropy X\\/x± = 10 ^- The symmetric difference scheme oF is founded 
on the support operator method (SOM) introduced in Reference^®, which is also known as 
Mimetic Finite Difference (MFD) method. Its main idea is to dehne the divergence operator 
as the adjoint of the gradient operator, resulting in a self-adjoint diffusion operator. This 
concept has also reached LDG methods, which are effectively used to solve parabolic and 
elliptic equations^^. However, so far the properties of LDG method have not been investi¬ 
gated for anisotropic heat conduction in magnetized fusion plasmas. 

Aligned schemes rely either on or 3D interpolation^®, and fix either the toroidal 

angle /S.ip or the length of the held line As. The interpolation points are found by tracing 
the magnetic held line with an appropriate integration method. Here, less accurate schemes 
neglect the detailed curvature of the magnetic held line and employ an approximate aligned 
stencil. 

The hux-coordinate independent (FGI) approach is an aligned 2D interpolation scheme, 
which approximates the parallel gradient Vy with a central hnite diherence®^. In that case, 
high-order 2D interpolation techniques have ehectively reduced numerical perpendicular dif¬ 
fusion. Here, a parallel advection equation supported the usability of the FGI approach even 
for diverted axisymmetric magnetic helds®^. The FGI approach was adapted to the MFD 
method to tackle the anisotropic dihusion equation®. It was shown that the self-adjoint dis¬ 
cretization of the dihusion operator exhibits less numerical perpendicular dihusion than a 
naive discretization®’^’®®. In this paper we lay the focus on the numerical investigation of the 
convergence rate of the aligned hnite diherence schemes for magnetic helds with V ■ h 7 ^ 0, 
which has not yet been considered. 

Aligned diherences with 3D interpolation can be extended to yield comparable numerical 
perpendicular heat huxes to non-aligned hnite diherences in the case of magnetic helds with 
V ■ b = 0®®. However, for V ■ 6 7 ^ 0 and closed held lines along the symmetry axis it was 
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unfolded that convergence is lost for high levels of anisotropy. 

The emphasis of all the so far mentioned schemes is put on magnetized fusion plasmas where 
only low temperature gradients along the magnetic held are expected. An extensive discus¬ 
sion on monotonicity preserving schemes, which are eligible for astrophysical plasmas can 
be found in Reference^®. 

In this work we present a non-aligned scheme based on the local discontinuous Galerkin 
(LDG) method. This scheme is formulated in the spirit of the SOM^®, so that the resulting 
diffusion operator is self-adjoint. This results in exact energy conservation and a vast redac¬ 
tion of unphysical perpendicular heat huxes. Moreover the order of the scheme is arbitrary. 
However, the non-aligned discretization suffers from a restrictive Gourant-Friedrichs-Lewy 
(GFL) condition, which is circumvented by an implicit time integrator. On the other side we 
implemented two aligned schemes of second order into the dG method, which were motivated 
from hnite differences and rely on 2D interpolation^’^^. The hrst scheme is not self-adjoint 
and fails in conserving energy exactly. We show nnmerically that in contrast to the direct 
aligned FD scheme the convergence is lost for higher resolutions. Hence the direct aligned dG 
scheme is not robust for general axisymmetric magnetic fields. The second aligned scheme 
is also self-adjoint so that energy is conserved. This scheme is corrected dne to jnmps at the 
cell bonndaries, which are accompanied by an unfavourable GFL condition and require an 
implicit treatment. Nnmerical tests reveal irregnlar and rednced convergence rates, which 
is in accordance with the resnlts of the self-adjoint aligned FD discretization. This is owed 
to corrugations which occnr on the grid scale. However, both aligned dG schemes demon¬ 
strate low perpendicular pollution even in cases of coarse resolntions and are particularly 
advantageous for complex flute-modes with k\\ ^ 0. All our proposed schemes are tested for 
an axisymmetric magnetic held where the divergence of h is not vanishing. The compnta- 
tions of the dG schemes were conducted with FELTOR (Full-F ELectromagnetic model 
in TORoidal geometry), which exploits the high degree of parallelism on GPUs and GPUs 
on shared and distribnted memory systems^^. GRILLIX provides the simnlations for the 
aligned FD schemes via an efficient hybrid parallelization for GPUs^®. 


4 



II. THEORY 


The theory part starts with an introduction to the LDG method in one dimension in 
Section II A. The generalization to three dimensions in cylindrical coordinates can be used 
for the non-aligned discretization of the parallel derivative in Section IIB. In Section IIC we 
present the aligned discretizations of the parallel derivative. Those make use of the Legendre 
polynomials for the 2D interpolation. The theory part is closed with a short description of 
the implemented time-discretizations in Section IID. 

A. The LDG method 

In the LDG framework one rewrites the model elliptic equation 

as a set of hrst order partial differential equations 


9 = 

dx 

in 

Q 

(5a) 

9 = 

= X9 

in 

Q 

(5b) 

dq 

dx 

= P 

in 

Q 

(5c) 


with a given function x{^) > 0 on D and corresponding boundary conditions on dQ, which 
are specihed later. In an LDG discretized scheme this set of equations decouples^ . After a 
short introduction to Legendre polynomials in Section IIA 1 we will derive the discretization 
of the hrst Eq. (5a) in Section IIA 2. With the help of the adjoint of the hrst derivative 
in Section IIA 3 we obtain a discretization of Eq. (4) in Section IIA 4. In the following we 
will retain the discussion to one dimension since the generalization to more dimensions is 
straightforward on orthogonal grids. 

1. Legendre polynomials 

Our discontinuous Galerkin methods are based on the use of orthogonal Legendre poly¬ 
nomials up to order P — 1. We dehne and wj, j = 0,..., P — 1 denoting the abscissas and 
weights of the Gauss-Legendre quadrature on the interval [—1,1] and denote Pk{x) as the 


5 



k — th Legendre polynomial (see e.gd®). With this, we define the forward transformation 
matrix 




( 6 ) 


Let us now consider an interval [a, b] and an equidistant discretization by N cells with 
cell center Xn and grid size h = in addition, we set := Xn + Given a function 
/ : [a, 6] —)■ M we then define fnj ■= and denote its dG expansion 

N P-1 

= (7) 

n=l fc=0 

with := J2j=o fnj- The expansion Eq. (7) is discontinuous across cell boundaries, 
hence the name dG. Here, Pnk{x) is the k — th Legendre polynomial in cell n 


Pnk{x) 


{ Pk {1{X - Xn)) , for X - Xn e |] 

0 , else. 


Let us define the diagonal weights matrix W and its inverse V 


( 8 ) 


H- := 

= if-/'’- 


(9a) 

(9b) 


The use of Legendre polynomials yields a natural approximation of the scalar product 
via Gauss-Legendre quadrature 

N P-1 


{fh-i 9h) ■ 


II/."" 


h\\L^ 


:= [ fhQhdx = EZ ^InjOnj =■ f^(l G> fT)g, 

n=l j=0 

w P-1 , 

n=l 7=0 


(10a) 

(10b) 


where 0 denotes the Kronecker product and the vector f contains the coefficients fnj- With 
these formulas we have a simple, accurate, and fast method to evaluate integrals. This 
is applied, for example, to compute errors in the L^-norm. Function products are easily 
computed pointwise, i.e. 


ifhQh) ni fnidn 


( 11 ) 
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2. The first derivatives 


The derivative of Pni{x) is not well defined on cell boundaries, which is why we now retain 
to a weak formulation. Consider 


dxfh{x)pni{x) dx = fhP, 


'Cn 


^n + l /2 


'Cn 


fh{x)dxPni{x) dx 


( 12 ) 


on the bounded cell domain Cn = \xn-i/ 2 -iXn+i/‘ 2 \ on 

The approximation fh{x) is double valued on the cell boundaries. The boundary terms have 
to be replaced by 


N P-1 


9 h{x)pni{x) dx = fp, 


'Cn 


I ^n + 1 /2 

m|a^n-l/2 


r^PnkdxPniix) dx, 

;i,^Q 


(13) 


where f{x) is the numerical flux across cell boundaries and we call gh{x) the numerical 
approximation to the hrst derivative. Let us dehne 


/+(x) := \im fh{x + e), (14a) 

£-5>0,£>0 

f~{x):= \im fh{x-e). (14b) 

e^>0,e>0 

We construct the numerical flux by the mean value {{/}} := (/'*' + /“) /2 and the jumps 
[[/]] := (/■•■ — /“) /2 at the cell boundaries 


/ = {{/}}+C||/]| 


(15) 


with C := {0,1, —1} determining the centered(C), forward(F) and the backward(B) flux 
respectively^®. For /: [a, 6] —)■ M we assume that 


f{h +e) = f{a +e), f{a-e) = f{b-e), (16a) 

/(a) = f{b) = 0, (16b) 

/(a) = /^(a)i f{b) = f~{b). (16c) 


for periodic, homogeneous Dirichlet and homogeneous Neumann boundary conditions re¬ 
spectively. Depending on what flux we choose, we arrive at various approximations to the 
derivative, e.g. for P = 1 (i.e. a piecewise constant approximation in each cell) our scheme 
reduces to the classic centered, forward and backward hnite difference schemes respectively. 
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The derivative can be written as a matrix-vector multiplication 


g = DJ. (17) 

For a concise notation we define 

:= (1 0 1/) ■ (1 0 F'^) ■ ■ (1 O F). 

We now show the matrix representation of the centered, forward, and backward one¬ 
dimensional discrete derivative for periodic boundary conditions that can be used in the 
implementation 

(M - M^) RL -LR \ 

-LR RL 

-LR ... , (18) 

RL 

RL -LR 





which we introduce mainly for ease of implementation. If a block-matrix class is written 
and the operations -t-, — and * are defined on it, the assembly of the derivative matrices is 
simplihed to a large extent. Note that for P = 1 we recover the familiar finite difference 
approximations of the first derivative. Finally, we note the boundary terms for homogeneous 

TABLE I: Upper left and lower right matrix entries for various boundary conditions. For 
Dirichlet and von Neumann BC the upper right and lower left entries are zero. 



DJ (forward) 

Dj. (backward) 

(centered) 


left 

right 

left 

right 

left 

right 

periodic 

-(M + L)T 

-(M + L)T 

(M + L) 

{M + L) 


\{M -M'^) 

Dirichlet 

-AfT 

-(M + L)T 

(M + L) 

-MT 

|(M-MT + L) 

\{M -M'^ -R) 

Neumann 

-(M + L)T 

M 

M 

{M + L) 

\{M -M'^ -L) 

\{M - M'^ + R) 


Dirichlet and Neumann boundaries in Table (I) noticing that only the corner entries of the 
matrices change. 

3. The adjoint of a matrix 

Recall that the adjoint of a matrix A is defined by the scalar product, i.e. 

■ [(1 ® IF) ■ A] ■ g = gT ■ ■ (1 ® IF)] ■ f =: g^ ■ [(1 ® W)A^] ■ f. 

From here we immediately get the relation 

At = (l(g)F)-A^-(l(g)IF). (22) 

There is a close connection between symmetric, A = A"*", and self-adjoint, A = Af matrices. 
If and only if the matrix A is symmetric, then (1 ® F) ■ A is self-adjoint. Of course, we have 
(A ■ P)t = Pt ■ At and (At)t = A. 

4 . Discretization of elliptic equations 

We discretize the one-dimensional general elliptic equation Eq. (4) on the interval O := 
[a, b]. We either choose periodic, Dirichlet, or Neumann boundary conditions on the left and 
right border for /. Our plan is to simply use one of the discretizations developed in the 
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last section for the right derivative and its adjoint for the left derivative to obtain an overall 
selfadjoint scheme. The multiplication with x is a simple function product 

Dl-x-DJ = p, (23) 

where is either D^, D~, or with the correct boundary terms and Dl is its adjoint 
according to Eq. (22). Note, that^® originally only proposed to use the forward or backward 
discretization for D^- Eq. (23) is indeed a self-adjoint discretization for the second derivative. 
However, it turns out that Eq. (23) is inconsistent for given p{x). The solution does not 
converge for P > 1. This problem is solved according to^® by adding a jump term [[/]] to 
the numerical flux of q in the case P > 1, which enhances the stability and accuracy of the 
LDG method 


? = {{9}}-c Ml-[[/]]. 


(24) 


That means we have to alter our discretization (23) according to 


where again J 


{1(^V){10 F) and 

/ {LL + RR) -RL 


J = 


-LR 


(LL + RR) -RL 
-LR 


-LR \ 


-RL 


(25) 


(26) 


\ -RL -LR {LL + RR)J 

for periodic boundaries. Again we give the correct boundary terms for Dirichlet and Neu¬ 
mann boundary conditions in Table II. Note that J is symmetric, thus the overall dis¬ 
cretization remains self-adjoint. Indeed, with Eq. (25) recovers the discretization 

proposed by^®. In addition, we remark that the centered discretization in Eq. (25) is symmet¬ 
ric with respect to an inversion of the coordinate system x —)■ —x even for double Dirichlet 
or Neumann boundaries, while the forward and backward discretization is not. 


B. Non-aligned discretization of the parallel derivative 

The generalization of the LDG method to higher dimensions for orthogonal grids is 
straightforward. All the matrices derived above can readily be extended via the appro- 
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J 


left 

right 

periodic 

L + R 

L + R 

Dirichlet 

L + R 

L + R 

Neumann 

R 

L 


TABLE II: Top left and bottom right entries for jump matrix 


priate Kronecker products. The space complexity of the matrices derived is 0{P‘^N) in one, 
0{P^N‘^) in two and 0{P'^N^) in three dimensions. In our case we use cylindrical coordi¬ 
nates {i?, Z, (p}, so that the volume form reads RdRdZdt^. All our presented dG methods 
are discretized on a cylindrical mesh, which consistis of NrNzN^ cells with PrPzOlp poly¬ 
nomial coefficients. This has to be taken into account by multiplying the weights in Eq. (22) 
by the factor R when computing the adjoint of matrices. We do so for the parallel deriva¬ 
tive, which can be written as Vy = + h^dz + The discrete centered, forward and 

backward derivatives are 


b^oD% + b^ oD^z + b"^ ° Dl: 

(27) 

b^oD+ + b^oD+ + b'^o D+, 

(28) 

b^oD^ + b^ oDz+b^oD-, 

(29) 


where o denotes pointwise multiplication. 


C. Aligned discretization of the parallel derivative 


We begin with the formulation of a held-ahgned discretization. If s denotes the held line 
following coordinate, then the one-dimensional discrete derivative along the held line reads 

^ J'k-\-l J'k —1 (30) 

From diherential geometry we know that to every smooth vector held b there is a unique 
curve of which the tangent in a point is the value of b at that point^'^. It is given by the 
solution of the diherential equation 


dz^ 

ds 



(31) 
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where is one of (i?, Z, (p) and 6 * are the contravariant components of h in cylindrical 
coordinates. Moreover, by dehnition we have 


d/(z(g)) 

ds 


= b ■ V/U,., 


(32) 


along a held line parameterized by distance s, i.e. instead of Vii/ we can choose to discretize 


ds * 


Let us divide the (p direction into equidistant planes of distance Acp. Unfortunately, 
from Eq. (31) we cannot easily determine the distance As for given Acp. It is better to 
integrate 


dz* 6* 5* 

dt ~ ~ 

since in this case d(p/dt = 1 ^ t = (p. We get 

dR _ 
dtp B'fi ’ 
dZ _ B^ 
dp) B'^ ’ 

together with the equation 

ds 1 i? 
dp 

for the length of the held line s. Eqs. (34) are integrated from p = d to p 
characterize the how generated by h/h'^ by 


(33) 


(34a) 

(34b) 


(34c) 

dzAp. We 


;= Z, p] := {R{±Ap), Z{±Ap), p ± Ap), (35) 

where {R{p), Z{p), s{p)) is the solution to Eqs. (34) with initial condition 


(/2(0),Z(0),s(0)) = (i?,Z,0). (36) 

Obviously we have = 1, but is not unitary since is not divergence free. 

The proposed centered discretization (30) for the parallel derivative then reads 

„ . ^ ^ ^ . / dpd - f ,„ 7 , 

'' ds dip ds s(+A^) — s(—Aip) ’ 

which is slightly diherent from Reference^^, where the relation (34c) was used to replace 
dp/ds. Previous derivations needed to construct an interpolation scheme in order to eval¬ 
uate functions on the transformed coordinates, which in general do not coincide with grid 
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points. Note that / is given as a dG expansion in i?, Z, (p with Pr,Pz arbitrary and P^p = 1. 
Hence the interpolation of the transformed points is naturally given by Eq. (7) (the 

extension to three dimensions is immediate). If {Rni, Zmj,p>k) are the grid points, we call 
:= TXp[Rni,Z^j,p>k] and (i?“ , := <^fc] the trans¬ 

formed coordinates along the held lines. We then have 

Nr Nz Pr—1 Pz —1 

/(Tiz) = /(fl+,z+j,^,+.)= x; x; E E 

mR=lmz=l Ir= 0 lz=0 

N<p Nr Nz Pr — 1 Pz — 1 

\^ (j+\'m,^mRjRmzjz r 

/ , / , / , / , G Jknimj J m^rnRjRmzjz ^ 

m^=l mR=l mz=l jR=0 jz=0 

(38a) 

Nr Nz Pr—I Pz — ^ 

/(Ti„z) = /(fl-,Z-j,^»_,)= E E E 

mR=lmz=l Ir= 0 lz=0 

Nip Nr Nz Pr — 1 Pz — 1 

( T-\m^rriRjRmzjz f 

/ , / , / , / , G tknimj J m^mRjRmzjz i 

m^=l mR=l mz=l jR=0 jz=0 

(38b) 


where the backward transformations f := (1 ® F)i are hidden in /. Thus, the interpolation 
of all the necessary points can simply be written as a matrix-vector product, where the 
interpolation matrices and I~ are independent of time since the magnetic held is assumed 
constant in time. The order of this interpolation is given by P, the number of polynomial 
coefficients. A consistency check is the relation ■ /“ = 1. 

The centered discretization (37) can now be written as a matrix vector product according 
to 


Qof = So ■ [/+ - /-] /, (39) 

where Sq is the diagonal matrix that contains the entries (^(-l-Ay)) — s(—A(/9))~^. This dis¬ 
cretization is not skew-symmetric since the held lines are not volume-preserving, or 7 ^ 

I~. The forward (-I-) and backward (—) operators Q for the parallel derivative are formulated 
analogously to 


Q+f 

Q-f 


s+ ■ [/+ - /] /. (40a) 

S- ■ [/- - /] /, (40b) 


where S± is the diagonal matrix that contains the entries (^(iA^?)) 


1 
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1. Boundary conditions 


The main problem with the above scheme is the question what to do when a held line 
crosses the simulation boundaries. Boundary conditions(BC) are formulated in cylindrical 
coordinates and are either homogeneous Dirichlet or Neumann. The latter are straightfor¬ 
ward and implemented by setting Vy / = 0 if the held line leaves the computational domain. 
For Dirichlet BCs one idea is to simply cut the contribution from held lines that leave the 
computational domain. While this works in practice it is unclear what numerical and physi¬ 
cal side-ehects this procedure might have. Another idea is to check to every point z whether 
Ta<^z lies inside our simulation box or not. If not, we have to hnd where exactly the held 
line intersects the simulation box. We have to hnd ifb such that the result of the integration 
of Eq. (34) from 0 to lies on the boundary. The angle (pb can be found by a bisection 
algorithm knowing that 0 < (pb < This kind of procedure is known as a shooting 

method. When all points are found, ghost cells can be constructed in the correct way. 

D. Time-stepping 

Having discretized the spatial derivatives of Eq. (3), we can now independently choose 
a discretization in time. This is known as the method of lines. In order to be able to test 
both explicit and implicit methods we implemented a semi-implicit Runge-Kutta (SIRK) 
scheme^^. The scheme is of third order in the explicit terms and of second order in the 
implicit terms. The method combines the cost-effectiveness of explicit algorithms with the 
stability of purely implicit time-steppers. We solve the implicit substeps by a conjugate 
gradient method, which works as long as the implicit part remains symmetric and linear. 
We hnally remark that, given the explicit and implicit parts, an implementation of the time- 
stepper only require linear algebra vector space operations. In fact, all of the up to now 
discussed algorithms can be implemented using basic vector-vector operations and a sparse 
matrix-vector product. 

III. THREE DG SCHEMES 

We are now ready to propose three different numerical approximations D to the parallel 
Diffusion operator An. 
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A. Direct aligned dG discretization 


In the first approximation we apply the divergence free property of the magnetic held B 
and write V ■ (bf) = -SVy We then use one of the centered, forward or backward 

aligned discretizations of the parallel derivative of Section IIC to get a discretization of 
A||/ = 5V|| (5“iV|j/). Hence, we can formulate two schemes based on either centered 
differences or an arithmetic average over the forward and backward differences 

D^f = BQoB-^Qof, D^f = ^ [BQ+B-^Q+f + BQ_B-^Q_f] . (41) 

In contrast to the naive discretization of Reference^, this scheme holds for magnetic helds 
with V ■ b ^ 0 and makes use of the Legendre polynomials for the interpolation. 

B. Self-adjoint aligned dG discretization 

We note the adjoint of the parallel derivative: 

Vj/ = -V ■ (i /) . (42) 

With this relation we can dehne the parallel diffusion operator of Eq. (2) to be self-adjoint 
according to 

A||=-VjV||. (43) 

This leads to the idea to simply use a discretization for the parallel derivative and use the 
matrix adjoint Eq. (22) to get a discretization for the divergence. Hence, we can derive the 
self-adjoint discretization for centered and averaged differences 

= [qIQo + j] /, D^f = \ [QiQ+ + Q^-Q- + 2j] /. (44) 

The scheme is similar to that of Reference^ except for the interpolation technique, and that 
in the dG framework we have to add the jump term J of Eq. (25) to obtain a convergent 
and stable diffusion operator. 

C. Self-adjoint non-aligned LDG discretization 

The directional derivative Vy = 6 ■ V can be discretized directly as Vy = h^dn + b^dz + 
using the dG methods developed in the Section HA2. The adjoint of the resulting 
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matrix is then a discretization for the divergence V ■ {bf) again. The compostion of both 
yields then a discretization for the parallel Laplacian according to 

DSf = [gJG„ + j] /, D°J = t [gVG+ + GlG_ + 2j] /. (45) 

Note, that we add jump terms for stabilization once more. 

IV. TWO ALIGNED FD SCHEMES 

The dG schemes are compared to two aligned FD discretizations, presented in Refer- 
ences^’^’^®. The first aligned FD scheme is a generalization of the naive scheme, which is 
based on third order bipolynomial interpolation (N-3) and was originally proposed to dis¬ 
cretize Vy. The operator Vy was added to this scheme so that effects arising from 

V • 6 7 ^ 0 are taken into account. The discrete parallel gradient is thereby computed via 
a central hnite difference along magnetic held lines^^, and the resultant scheme is denom¬ 
inated D-3. The second aligned FD schemes S-3 is formulated in the spirit of the SOM. 
In the S-3 scheme the discrete parallel gradient is also computed with the help of third 
order bipolynomial interpolation. The parallel diffusion operator then follows via requiring 
its self-adjointness property on the discrete level. 

V. TESTS 

We use an axisymmetric magnetic held with V-h 7 ^ 0. In cylindrical coordinates {R, Z, ip} 
it is given by 

B = RVip + V'lpp X Vp>, (46) 

with Jo = 20 , Rq = 10 and a poloidal hux of the form 

'ipp = cos {^{R - Ro)^ cos . (47) 

We perform all our numerical tests on the domain R G [9,11], Z G [—1,1] and ip G [0, 27r], 
except for the test given in Section VB 3. With this choice of parameters the magnetic held 
never crosses the simulation box. 
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A. Convergence test to an analytical test function 


First, we define a general test function / = —ipplR, Z) cos {(p), which fulhlls Dirichlet 
as well as Neumann boundary conditions. We use this function to directly compute the 
relative error eAj|/ = \\Df — A||/||i 2 /||A||/|| 2,2 between the numerical diffusion operator 
Df and the analytical solution A||/. At this point we introduce also a second relative 
error e/ = ||/b - /HlVII/IU^, which is computed by inverting the diffusion operator with 
a conjugate gradient method. On the continuous level this would result in an ill-posed 
problem. However, the discrete parallel diffusion operator contains always a small amount 
of numerical perpendicular diffusion, which makes the discrete operator invertible. 

All the so far mentioned aligned schemes are constructed to converge with order two in 
Aif for a homogeneous magnetic held^. Moreover, the order of the non-aligned scheme is 
arbitrary. We dehne this as the optimal order. For inhomogeneous magnetic helds we can 
observe this optimal order only if we increase the resolution in A/j, Nz and with the 
correct factors, which depend on the interpolation order and the order in (/^-direction. Our 
tested schemes are limited to Prz = 3 and = 1. For the non-aligned LDG method we 
additionally investigate the convergence rates for = 3 and P^ = 3. 


1 . dG schemes 

We are now ready to show the orders of convergence of the centered and averaged dis¬ 
cretizations. For this purpose we increase the resolutions in all directions with the correct 
order according to Prz = 3, Nr = Nz = 5 (2*^/^), N^ = 5 (2*) with i = {0,1,..., 5}. The 
relative errors and orders of convergence are depicted in Table III for the centered schemes 
and Table IV for the averaged schemes. Those tables reveal that the non-aligned LDG 
discretization is the only scheme, which has the optimal order of convergence for P^^ = 1 and 
also superconvergent rates for P^ = 3. The aligned schemes do both show up stagnating 
convergence rates for higher resolutions whereas the self-adjoint aligned scheme additionally 
fails in converging with the optimal order of two. Instead it converges with a reduced and 
irregular order, casting it into a true supraconvergent discretization. This is a consequence 
of oscillatory relative error helds, which are shown in Figure 1 and discussed throughout the 
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TABLE III; error and orders of accuracy 

(orderA</, := log(L^error^^-/L^errorA(/,)/log(A0“/A(;/))) for the centered non-aligned LDG 
scheme, the centered direct aligned scheme and the centered self-adjoint aligned scheme. 
The error for the self-adjoint schemes is computed with e/ and for the direct aligned 

scheme with eA,,/ 



= 1 


CO 

II 

A(/9 

error order 

nD 

error order 

n-4 

L?' error order 

A<^ 

error order 

27r/3 

4.85E-F00 - 

8.29E-01 - 

6.65E-h00 - 

27r/3 

3.88E-02 - 

27r/5 

7.46E-01 3.66 

4.28E-01 1.30 

l.llE-hOO 3.51 

27r/5 

7.58E-03 3.20 

27r/10 

1.43E-01 2.39 

1.25E-01 1.77 

3.39E-01 1.71 

27r/10 

7.95E-04 3.25 

27r/20 

3.36E-02 2.09 

3.26E-02 1.94 

1.92E-01 0.82 

27r/20 

6.77E-05 3.55 

27r/40 

8.27E-03 2.02 

8.24E-03 1.98 

1.25E-01 0.62 

27r/30 

1.61E-05 3.55 

27r/80 

2.08E-03 1.99 

2.09E-03 1.98 

4.50E-02 1.48 

~ 

- 

27r/160 

5.14E-04 2.02 

5.94E-04 1.82 

6.00E-02 -0.42 

- 

- 


next sections. Furthermore, we note that we observed stagnating convergence rates also for 
the relative error of the parallel gradient operator (5{o,+,-}- 


2. Aligned FD schemes 

Additionally the same convergence test was performed with the GRILLIX code, which 
is based on the aligned FD schemes depicted in Section IV. 

In Table V the relative error eAy/ and convergence rates for the aligned FD schemes are 
shown. We adapted the resolution according to roughly Nji = Nz = 3 ■ 5(2*^/^), = 

5(2*), i = {0,1,... 6}, where the additional factor of 3 at Nji, Nz shall compensate for the 
factor Prz = 3 used in the test of Section VAl. The derived relative error for the D-3 
scheme is of order two. However, similar to the scheme it does not conserve the self¬ 
adjointness property of the parallel diffusion operator and therefore does not conserve the 
norm, i.e. energy (see also Section VBl). The S-3 scheme exhibits an irregular and 
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TABLE IV: error and orders of accuracy 

(order:= log(L^error^^-/L^errorA</,)/log(A0“/A(;/))) for the averaged non-aligned LDG 
scheme, the averaged direct aligned scheme and the averaged aligned self-adjoint scheme. 
The error for the self-adjoint schemes is computed with e/ and for the direct aligned 

scheme with eA,,/ 



= 1 


CO 

II 


■^a 




■^a 

Acp 

L?' error order 

error order 

error order 

Acp 

error order 

27r/3 

4.62E-01 - 

3.17E-01 - 

7.77E-01 - 

2n/3 

3.95E-02 - 

2tt/5 

1.43E-01 2.30 

1.25E-01 1.82 

2.13E-01 2.53 

27r/5 

7.96E-03 3.14 

27r/10 

3.36E-02 2.09 

3.26E-02 1.94 

6.66E-02 1.68 

27r/10 

8.46E-04 3.23 

27r/20 

8.27E-03 2.02 

8.24E-03 1.98 

4.43E-02 0.59 

27r/20 

8.94E-05 3.24 

27r/40 

2.06E-03 2.01 

2.17E-03 1.92 

2.42E-02 0.87 

27r/30 

2.42E-05 3.22 

0 

00 

CNl 

5.15E-04 2.00 

7.69E-04 1.50 

2.29E-02 0.08 

- 

- 

27r/160 

1.29E-04 2.00 

5.05E-04 0.61 

4.46E-02 -0.96 

- 

- 


slower convergence behaviour. This is owed to corrugations (see Figure 1), which occur if 
the transformed coordinates [i?j, Zj, (p^] of two neighbouring grid points are in the same 
grid cell, respectively jump across grid cells. The mechanism is described in more detail in 
Reference^. Apart from these corrugations the error is of the same level as for the D-3 
scheme. In general, also for magnetic helds with V ■ 6 = 0 corrugations can occur and have 
been observed, but they are more pronounced in situations with V ■ h 7 ^ 0. As long as the 
corrugations remain on the grid scale (a suitable criterion for this is given in Reference^) 
they could possibly cured by adding a small amount of perpendicular diffusion, which is 
usually anyway present in turbulence simulations. 


B. Numerical simulations 

We investigate now the conservation properties and the convergence behavior of our 
presented schemes in numerical simulations of Eq. (3). In this respect this we use Neumann 
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TABLE V: Relative error eAy/ and its order of accuracy 
(orderA <)4 := log(L^error^^-/L^errorA</,)/log(A(;/)“/A(;/))) for the aligned FD schemes. 


A(p 

Nr^z {h) 

D-3 

error order 

S-3 

error order 

2ttI3 

20 (l.OOE-01) 

3.06E-01 - 

3.06E-01 - 

I'Kjh 

20 (l.OOE-01) 

1.20E-01 1.81 

1.42E-01 1.49 

27r/10 

24 (8.33E-02) 

3.28E-02 1.89 

1.16E-01 0.30 

27r/20 

40 (5.00E-02) 

5.12E-03 2.00 

4.42E-02 1.39 

27r/40 

60 (3.33E-02) 

2.12E-03 1.95 

2.49E-02 0.83 

27r/80 

100 (2.00E-02) 

5.12E-04 2.04 

1.65E-02 0.59 

27r/160 

150 (1.33E-02) 

1.25E-04 2.04 

1.24E-02 0.41 

27r/320 

240 (8.33E-03) 

2.81E-05 2.14 

9.35E-03 0.41 



FIG. 1: The error helds Df — Ay / of the convergence test are shown for the self-adjoint 
aligned dG scheme and the self-adjoint aligned FD scheme S-3 at = 0. The 

resolution is = 1, Prz = 3, = 80, Nr = Nz = 32 and equally 

= 80, Nr = Nz = 100 for the latter. In contrast to the S-3 scheme the 
oscillations/corrugations are more pronounced for the scheme. 


boundary conditions to show explicitly the proposed energy conservation of our self-adjoint 
schemes. The convergence is examined for a general non-aligned three dimensional Gaussian 
blob in Section VBl. The maintenance of the most basic k\\ = 0 mode is shown for all 
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schemes in Section VB. We take X|| = 100 and = 0 for all of our simulations. The 
time integration was either performed with a second-order implicit or third-order explicit 
Runge-Kutta scheme. The time step was fixed to At = 0.001 so that the time error can be 
ignored. The convergence test was performed with an even smaller timestep. 


1. Gaussian Blob 


We initialize a non-aligned Gaussian blob in ip direction of the form 

{R-RhY {Z-Z,)^ 


fo = A exp 


24 


24 


24 


(48) 


with A = 0.1, Rh = 10.6, Zf, = 0.0, pi, = tt, = az = 0.1, and = 0.5. To show the 
convergence of our schemes in time simulations for a general three dimensional perturbation, 
we perform a short highly accurate simulation over 100 time steps to t = 0.01 with the non- 
aligned LDG scheme for a resolution of Prz = 3, Nr = Nz = 51, = 1 and = 162. 

This numerical reference state /i is used for the computation of the relative error = 
Wft - /llUa/il/lll i 2 at t = 0.01. The orders of convergence are listed in Table VI. The high 
perpendicular resolution should lead to convergence rates of order two if only the resolution 
in p is incremented. We obtain second order convergence for the non-aligned LDG scheme 
and also for the direct aligned discretization. However, the self-adjoint aligned discretization 
has again reduced orders of convergence. This is contributed to the oscillatory temperature 
helds, which are depicted in Figure 2 for a blob with A = 0.1, Rh = 10.6, Z^ = 0.0, pi, = n, 
= o'z = 0.1, and = 2. The oscillations arise because the adjoint of the interpolation 
matrices do no longer coincide with the three dimensional weights for all cells. Hence, 

the transformation is not unitary since V ■ h 7 ^ 0 for the investigated magnetic field. For 
the aligned FD schemes the position of the the interpolation points is always in the center 
of the sampling points of the interpolating polynomials. As a consequence the corrugations 
are better suppressed than for the aligned dG schemes. The initialized structure is not held 
aligned and bean like temperature helds appear for all aligned schemes, whereas the LDG 
scheme yields smooth results, which is depicted in Figure 3. 

In Figure 4a) and Figure 4b) we show that energy is conserved for all self-adjoint schemes. 
The direct aligned schemes fail in conserving the norm. We note here that the self-adjoint 
schemes are also superior in terms of numerical stability. This is based on the fact that the 
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FIG. 2: These contour plots show the evolution of a three dimensional Gaussian blob at 
t = 1 and (p = TT for the non-aligned LDG scheme (left), the direct aligned scheme (center) 
and the self-adjoint aligned scheme (right). The latter shows up spurious oscillatory 
temperature helds. The bean like temperature helds of the aligned discretizations are a 
consequence of the hnite resolution in the direction of (p unlike the non-aligned 
discretization, which has smooth helds even for moderately resolved cases. The resolution 


was set to Prz = 3, Nr = Nz = 20, = 21. 

D-3 S-3 



9.5 10.0 10.5 11.0 9.5 10,0 10.5 110 

R R 


FIG. 3: Result of simulation of Gaussian blob at t = 1 and p = tt for aligned FD schemes. 

Resolution was set to Nr = Nz = 60, N^ = 20. 


2. Profile maintenance 

We initialize a prohle of the form /o = O.l'^p for which the parallel Laplacian vanishes 
exactly to zero Am/q = 0. Hence, the initial state is conserved over time and deviations 
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FIG. 4: a) The energy conservation over time for the averaged non-aligned LDG scheme 
(blue circles), the averaged direct aligned scheme (green rhombi) and the averaged 
self-adjoint aligned scheme (red squares). The resolution was set to Prz = 3, 

Nr = Nz = 20, N^p = 20. b) The energy conservation over time for the D-3 scheme (green 
rhombi) and the S-3 scheme (red squares). The resolution was set to Nr = Nz = 60, 

= 20 


TABLE VI; error and orders of accuracy 
(orderA,/, := log(L^errorA(^-/T^errorA<^)/log(A0“/A0)) for the averaged non-aligned LDG 
scheme, the averaged direct aligned scheme and the averaged self-adjoint aligned scheme at 

t = 0.01. The time step is hxed to At = 0.0001. 


Aip 

error order 

D? 

L? error order 

error order 

27r/6 

2.58E-04 - 

7.48E-04 - 

7.48E-04 - 

27r/18 

3.50E-05 1.82 

4.62E-04 0.44 

4.64E-04 0.44 

27r/54 

3.80E-06 2.02 

5.18E-05 1.99 

1.15E-04 1.27 

27r/162 


5.32E-06 2.07 

4.45E-05 0.86 


are a measure for the numerical perpendicular heat flux or numerical sinks and sources. 
Gonsequently, we compute the relative error via e/p = H/i — /o||l 2 /||/o||l 2 . In Figure 5a) 
we analyze its time evolution for three different perpendicular resolutions Nr = Nz = 
{5, 20, 80} hxing the resolution in (p to N^ = 10 and polynomial coefficients to = 1 
and Prz = 3. We observe that the aligned dG schemes are superior to the non-aligned 
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LDG scheme for low resolutions. However, for finer resolutions all three schemes do have 
nearly the same numerical perpendicular heat flux. This error decreases with the order 
of the interpolating polynomial, which is shown in Table VII at f = 5. Here, we observe 
again superconvergence for the non-aligned LDG scheme. Furthermore, we note that e/^ is 
independent of N^p for the non-aligned LDG scheme, whereas it increases with N'^ for the 
aligned schemes^. This is ascribed to the aligned differences, which are inversely proportional 
to for a fc|| = 0 mode. 

For the aligned FD schemes Figure 5b) shows that the D-3 scheme exhibits a much higher 
numerical perpendicular heat transport than the S-3 scheme, also if V ■ 6 7 ^ 0 (for details 
cf.^). 


a) b) 



FIG. 5: a) The relative error e/^ over time for the non-aligned LDG scheme (blue circles), 
the direct aligned scheme (green rhombi) and the self-adjoint aligned scheme (red squares) 
for three different perpendicular resolutions Nji = Nz = {5, 20, 80} (dashed, dotted, 
dot-dashed), b) The relative error ef^ over time for the D-3 scheme (green rhombi) and 
the S-3 scheme (red squares) for three different perpendicular resolutions 
Nr = Nz = {15,60,240} (dashed, dotted, dot-dashed). 


3. Aligned blob in tokamak geometry 

We verify our schemes now on a realistic axisymmetric tokamak equilibrium with lower 
X-point. The used Solov’ev equilibrium is an analytical solution of the Grad-Shafranov 
equation in cylindrical coordinates^^ and resembles to a sufficient extent a typical deuterium 
discharge in the GOMPASS tokamak with major radius Rq = 0.56m, inverse aspect ratio 
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TABLE VII; error e/p and orders of accuracy 
(orderAi? := log(L^errorA/j-/-L^errorAi?)/log(Ai?“/Ai?)) for the averaged non-aligned 
LDG scheme, the averaged direct aligned scheme and the averaged self-adjoint aligned 

scheme 


AR = AZ 

■^a 

error order 

^ a 

L?' error order 

error order 

2/5 

9.74E-04 - 

3.05E-04 ~ 

4.04E-04 - 

2/10 

7.22E-05 3.75 

3.11E-05 3.29 

6.22E-05 2.70 

2/20 

6.22E-06 3.54 

3.91E-06 2.99 

8.13E-06 2.94 

2/40 

6.44E-07 3.27 

4.89E-07 3.00 

l.OlE-06 3.00 

2/80 

7.31E-08 3.14 

7.75E-08 2.66 

1.36E-07 2.90 


e = 0.41, elongation k = 1.75, triangularity 6 = 0.47, toroidal magnetic held Bq = IT and 
reference electron Temperature T^o = 500eV^^. This hxes the drift scale pso = y/rriiTe/eBQ 
and the ion gyro-frequency flj = eBo/rrii, which normalize the spatial and time scales. The 
corresponding geometric coefficients are 


A = 0, 

C2 = -0.0866241743631724894540195833599, 
C4 = -0.0763123710053627236722106917166, 
C6 = -0.0915754123901870980044268057017, 
Cg = 0.0427189122507639110238025727235, 
cio = -0.130472413601776209642913874538, 
ci2 = 0.00421267189210391228482454852824. 


Cl = 0.0735011444550040364542354505071, 

C 3 = -0.146393154340110145566533995171, 

C 5 = 0.0903179011379421292953276540381, 

C 7 = -0.00389228297983755811356246513352, 
cg = 0.227554564600278024358552580518, 
cii = -0.0300697410847693691591393359449, 

(49) 


In Figure 6 we show the time evolution of an initially held aligned blob with /cy 7 ^ 0 at /? = 
3n/2 for a typical spatial resolution of Prz = 3, = 80, Nz = 120, = 1 and = 20. 

We initialize a two dimensional Gaussian blob at /? = tt of the form fo{R, Z, tt) (see Eq. (48)) 
and with parameters A = 0.1, Rb = i?o (1 + 0.8e), Zb = 0.0, ipb = n, = az = 5pso- With 
the help of the held line transformation 
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(t = 1000) 
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FIG. 6: These plots show the evolution of a held aligned blob from f = 0 to f = 1000 for 
the non-aligned LDG scheme (2nd column), the direct aligned scheme (3rd column) and 
the self-adjoint aligned scheme (4th column). Here, we show the toroidal cross section at 
ip = 37r/2. The non-aligned LDG discretization suffers from pollution by numerical heat 
fluxes, whereas the aligned dG discretizations are able to retain the initial blob structure 
whereby also highly sheared variants emerge in the left top corner of the numerical 
domain. The resolution was set to Prz = 3, = 80, Nz = 120, = 1 and = 20.. 

/o,|| (R, Z, TT ± A(p) = exp TllfoiR, Z, tt) (50) 

we transform the two dimensional intial held to all other planes. With cr^ = 0.57ri?o fhis 
yields a held aligned blob that is spread over approximately half of the toroidal domain. 
Due to the sheared magnetic held the two dimensional Gaussian blob structure is distorted 
at (p 7^ TT. As time progresses we expect that the initial cross section of the held aligned 
blob is conserved. Due to parallel heat conduction the held aligned blob traces out the 
magnetic held lines and highly sheared variants can emerge on the same toroidal plane. As 
depicted in Figure 6 this is not the case for the non-aligned LDG scheme. Here, the coarse 
resolution fails to resolve the complex flute-mode properly as time advances. The fc|| ^ 0 
flute-mode is dominated by numerical heat fluxes. Only the aligned dG schemes can resolve 
this mode at reasonably low resolution in p direction. The high numerical heat flux of the 
non-aligned LDG scheme is caused by the fact that held lines are sampled by the grid. In 
contrast aligned schemes resolve the magnetic held lines with a high accuracy based on a 
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highly accurate held line tracing procedure. We remark again the violation of the norm 
conservation of the direct aligned dG scheme in contrast to the self-adjoint schemes. For 
this choice of numerical parameters this violation becomes signihcant at t ~ 1800, which is 
shown in Figure 7. 



FIG. 7: The energy conservation over time of a held aligned Gaussian blob for the 
averaged non-aligned LDG scheme (blue circles), the averaged direct aligned scheme (green 
rhombi) and the averaged self-adjoint aligned scheme (red squares) is shown. 


VI. CONCLUSION 

We have presented three diherent dG discretizations for the parallel dihusion operator. 
The non-aligned LDG discretization is hexible in the order of convergence whilst keeping the 
numerical perpendicular heat huxes at a small level for a basic k\\ = 0 mode and conserving 
energy exactly. Our numerical measurements revealed that in general V ■ 6 7 ^ 0 cases 
the aligned dG schemes may exhibit stagnating, reduced or irregular convergence rates 
whereas the non-aligned LDG scheme remains robust. However, especially for hute-modes 
and coarse resolutions the aligned schemes are advantageous to the non-aligned LDG scheme. 
Most notably the non-aligned LDG discretization requires signihcantly higher resolutions for 
k\\ 7 ^ 0 flute-modes in typical X-point geometry. We compared the latter dG schemes to the 
aligned FD schemes of References'’^ for a magnetic held with V • 6 7 ^ 0. The direct aligned 
FD scheme converges with robust rates but is inferior to the self-adjoint aligned FD scheme 
regarding energy conservation and hute mode maintenance. On the downside the latter 
scheme showed up degraded and anomalous convergence rates due to the appearance of 
few corrugations on the grid scale. The herein presented dG discretizations for the parallel 
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diffusion operator are implemented into a Full-F ELectromagnetic model in TORoidal 
geometry (FELTOR) and support the computation of long term turbulence simulations of 
X-point equilibria. Those results are postponed to future publications. 
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